--- title: "3. Zonal statistics with discrete data" author: "Gabriel Carrasco & Antony Barja" date: 2020-12-02 tags: ["statisticzonal","discrete"] ---
The objective of this demo is to quantify the categories of herbaceous vegetation in San Diego city - USA. The dataset used is Copernicus Global Land Service (CGLS) which provides a series of bio-geophysical products on the status and evolution of land surface at global scale to a spatial resolution of 100m. More information click here
Information: |
- For this demo you need rgee, sf, tidyverse, mapview, raster and ggspatial packages previously installed. |
library(rgee)
library(sf)
library(tidyverse)
library(mapview)
library(raster)
library(ggspatial)
ee_Initialize()
> ee_Initialize()
── rgee 1.1.2 ────────────────────────── earthengine-api 0.1.292 ──
✓ user: not_defined
✓ Initializing Google Earth Engine: DONE!
✓ Earth Engine account: users/ambarja
───────────────────────────────────────────────────────────────────
sandiego <- st_read(
"https://github.com/healthinnovation/sdb-gpkg/raw/main/SanDiego_districts.gpkg"
) %>%
dplyr::select(PO_NAME)
head(sandiego)
> head(sandiego)
Simple feature collection with 6 features and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -117.1162 ymin: 32.58161 xmax: -116.208 ymax: 32.90116
Geodetic CRS: WGS 84
PO_NAME geom
1 Alpine MULTIPOLYGON (((-116.6075 3...
2 Bonita MULTIPOLYGON (((-116.9655 3...
3 Boulevard MULTIPOLYGON (((-116.2085 3...
4 Campo MULTIPOLYGON (((-116.3573 3...
5 Chula Vista MULTIPOLYGON (((-117.0728 3...
6 Chula Vista MULTIPOLYGON (((-117.0081 3...
dataset <- ee$Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019")$
select("discrete_classification")$
eq(30) # 30: Herbaceous vegetation
herbaceous <- dataset$updateMask(dataset)
ee_get_data <- function(img, region, scale = 1000) {
lista_histo <- list()
for (i in 1:nrow(region)) {
region_ee <- region[i, ] %>% sf_as_ee()
ee_histo <- img$reduceRegion(
reducer = ee$Reducer$frequencyHistogram(),
geometry = region_ee,
scale = scale
)
lista_histo[[i]] <- ee_histo$getInfo() %>%
map_df(., .f = as.data.frame) %>%
mutate(NAME = region[[i, 1]] %>% as.vector())
}
return(lista_histo)
}
ee_get_data# Time: ~ 2 min
rawdata <- ee_get_data(img = herbaceous, region = sandiego)
class(rawdata)
[1] "list"
# Herbaceous area in Hectare Units
newdata <- rawdata %>%
map_df(.f = as_tibble) %>%
mutate_if(
is.numeric,
.funs = function(x) {x * 1000 * 1000 / 10000}
) %>%
rename(area_Ha = X1)
head(newdata)
# A tibble: 6 × 2
area_Ha NAME
<dbl> <chr>
1 800 Alpine
2 1073. Bonita
3 188. Boulevard
4 1300 Campo
5 59.6 Chula Vista
6 179. Chula Vista
# Left join with spatial data
sandiego_herb <- left_join(
sandiego,
newdata,
by = c("PO_NAME"="NAME")
)
head(sandiego_herb)
Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -117.1162 ymin: 32.58161 xmax: -116.208 ymax: 32.90116
Geodetic CRS: WGS 84
PO_NAME area_Ha geom
1 Alpine 800.00000 MULTIPOLYGON (((-116.6075 3...
2 Bonita 1073.33333 MULTIPOLYGON (((-116.9655 3...
3 Boulevard 187.84314 MULTIPOLYGON (((-116.2085 3...
4 Campo 1300.00000 MULTIPOLYGON (((-116.3573 3...
5 Chula Vista 59.60784 MULTIPOLYGON (((-117.0728 3...
6 Chula Vista 179.21569 MULTIPOLYGON (((-117.0728 3...
mapview(sandiego,zcol="area_Ha", layer.name = "Herbaceous")